###################################################################################
# Replication files for: Negative References to Amicus Briefs in Judicial Reasoning
# Authors: Johan Lindholm, Daniel Naurin & Philipp Schroeder
# Published in: Journal of Law and Courts
# Contact: P.Schroeder@gsi.uni-muenchen.de

### File 3: This script provides the findings and figures for the robustness checks presented in the manuscript's supplementary materials

library(arm)
library(rstanarm)
library(broom)
library(tidyverse)
library(stringr)
library(pscl)
library(texreg)
library(sandwich)
library(lmtest)
library(ggplot2)
library(boot)
library(parallel)
library(margins)
library(lavaan)
library(lavaan.survey)
library(semTable)
library(nnet)

rm(list=ls())
homeFolder <- Sys.getenv("HOME")


# read data ---------------------------------------------------------------

load("data/mentions_issue_level.RData")
load("data/mentions_submission_level.RData")



# hurdle models for figure 5 in the supplementary materials ---------------

hurdle1 <- hurdle(mentioned_negative ~ cjeu_pos_auto.numerical + z.ms_pos_auto0 + z.ms_pos_auto1 + ag_cjeu_disagreement + ag_com_disagreement +
                    z.sources_case_law.ag + z.sources_primary.ag + z.sources_secondary.ag +
                    z.duration_since_ag + panel_size +
                    number_of_observations + legact_primary + legact_secondary +
                    administrative + agriculture + criminal + customs + environment + healthcare + intellectual_property +
                    labour + marketing + migration + procedure + procurement + no_area_of_law + social + tax + transport +
                    z.subject_matter.count + year_proceeding.factor,
                  data = mentions_issue_level, dist = "negbin")
summary(hurdle1)

coef.count1 <- summary(hurdle1)$coefficients$count
coef.zero1 <- summary(hurdle1)$coefficients$zero

hurdle2 <- hurdle(mentioned_negative ~ z.ms_pos_auto0 + z.ms_pos_auto1 + ag_com_disagreement +
                    z.sources_case_law.ag + z.sources_primary.ag + z.sources_secondary.ag +
                    z.duration_since_ag + panel_size +
                    number_of_observations + legact_primary + legact_secondary +
                    administrative + agriculture + criminal + customs + environment + healthcare + intellectual_property +
                    labour + marketing + migration + procedure + procurement + no_area_of_law + social + tax + transport +
                    z.subject_matter.count + year_proceeding.factor,
                  data = mentions_issue_level, dist = "negbin")
summary(hurdle2)

coef.count2 <- summary(hurdle2)$coefficients$count
coef.zero2 <- summary(hurdle2)$coefficients$zero

hurdle_point_estimates <- function(coefs, var_index){
  bin <- matrix(nrow = 3, ncol = nrow(coefs))
  bin[1,] <- coefs[,1]
  bin[2,] <- coefs[,1] - 1.96 * coefs[,2]
  bin[3,] <- coefs[,1] + 1.96 * coefs[,2]
  bin <- as.data.frame(bin)
  names(bin) <- row.names(coefs)
  
  # select coefficients
  out <- as.data.frame(matrix(nrow = 3, ncol = length(var_index)))
  names(out) <- names(bin)[names(bin) %in% var_index]
  out[1,] <- bin[1,names(bin) %in% var_index]
  out[2,] <- bin[2,names(bin) %in% var_index]
  out[3,] <- bin[3,names(bin) %in% var_index]
  
  out <- out[,var_index]
  return(out)
}

coef_plot.hurdle1.count <- hurdle_point_estimates(coef.count1, c("cjeu_pos_auto.numerical","z.ms_pos_auto0","z.ms_pos_auto1","ag_cjeu_disagreement1","ag_com_disagreement1",
                                                                 "z.sources_case_law.ag","z.sources_primary.ag","z.sources_secondary.ag","z.duration_since_ag",
                                                                 "panel_size","number_of_observations","legact_primary","legact_secondary"))

coef_plot.hurdle1.zero <- hurdle_point_estimates(coef.zero1, c("cjeu_pos_auto.numerical","z.ms_pos_auto0","z.ms_pos_auto1","ag_cjeu_disagreement1","ag_com_disagreement1",
                                                               "z.sources_case_law.ag","z.sources_primary.ag","z.sources_secondary.ag","z.duration_since_ag",
                                                               "panel_size","number_of_observations","legact_primary","legact_secondary"))



coef_plot.hurdle2.count <- hurdle_point_estimates(coef.count2, c("z.ms_pos_auto0","z.ms_pos_auto1","ag_com_disagreement1",
                                                                 "z.sources_case_law.ag","z.sources_primary.ag","z.sources_secondary.ag","z.duration_since_ag",
                                                                 "panel_size","number_of_observations","legact_primary","legact_secondary"))

coef_plot.hurdle2.zero <- hurdle_point_estimates(coef.zero2, c("z.ms_pos_auto0","z.ms_pos_auto1","ag_com_disagreement1",
                                                               "z.sources_case_law.ag","z.sources_primary.ag","z.sources_secondary.ag","z.duration_since_ag",
                                                               "panel_size","number_of_observations","legact_primary","legact_secondary"))


# plot coefficients for model 7 and 8
pdf(file = "figures/hurdle_issue.pdf", height = 12, width = 10)
par(mfrow=c(2,1))
# zero model
par(mar=c(5,18,3,3))
plot("", axes = FALSE,
     xlim = c(-1.75,1.75),
     ylim = c(0.75,13.25),
     xlab = "Logistic regression coefficient estimates with 95% confidence intervals",
     ylab = "", pch = 19, cex = 1.5,
     main = "Zero model: Negative reference (Models 7 and 8)", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-1.75,1.75,.25), cex.axis = 1.2)
axis(2,at = seq(13,1,-1), las = 1, tick = TRUE, cex.axis = 1.2,
     labels = c("CJEU position",
                "MS favouring no restrictions",
                "MS favouring restrictions",
                "Conflicting positions: AG/CJEU",
                "Conflicting positions: AG/Commission",
                "Sources cited by AG: Case law",
                "Sources cited by AG: Primary law",
                "Sources cited by AG: Secondary law",
                "Deliberation time",
                "Panel size",
                "Number of observations",
                "Issue concerns primary EU law",
                "Issue concerns secondary EU law"))
abline(h = seq(13,1,-1), lty = 2, lwd = .5, col = "grey")
# m1
points(coef_plot.hurdle1.zero[1,], 13:1 + .1, pch = 19)
segments(t(coef_plot.hurdle1.zero[2,]), 13:1 + .1,
         t(coef_plot.hurdle1.zero[3,]), 13:1 + .1, lwd = 2)
# m2
segments(t(coef_plot.hurdle2.zero[2,]), c(12,11,9:1) - .1,
         t(coef_plot.hurdle2.zero[3,]), c(12,11,9:1) - .1, lwd = 2)
points(coef_plot.hurdle2.zero[1,], c(12,11,9:1) - .1, pch = 24, bg = "white")
abline(v=0, lty = 2)
# count model
par(mar=c(5,18,3,3))
plot("", axes = FALSE,
     xlim = c(-1.75,1.75),
     ylim = c(0.75,13.25),
     xlab = "Negative binomial regression coefficient estimates with 95% confidence intervals",
     ylab = "", pch = 19, cex = 1.5,
     main = "Count model: Negative references (Models 7 and 8)", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-1.75,1.75,.25), cex.axis = 1.2)
axis(2,at = seq(13,1,-1), las = 1, tick = TRUE, cex.axis = 1.2,
     labels = c("CJEU position",
                "MS favouring no restrictions",
                "MS favouring restrictions",
                "Conflicting positions: AG/CJEU",
                "Conflicting positions: AG/Commission",
                "Sources cited by AG: Case law",
                "Sources cited by AG: Primary law",
                "Sources cited by AG: Secondary law",
                "Deliberation time",
                "Panel size",
                "Number of observations",
                "Issue concerns primary EU law",
                "Issue concerns secondary EU law"))
abline(h = seq(13,1,-1), lty = 2, lwd = .5, col = "grey")
# m1
points(coef_plot.hurdle1.count[1,], 13:1 + .1, pch = 19)
segments(t(coef_plot.hurdle1.count[2,]), 13:1 + .1,
         t(coef_plot.hurdle1.count[3,]), 13:1 + .1, lwd = 2)
# m2
segments(t(coef_plot.hurdle2.count[2,]), c(12,11,9:1) - .1,
         t(coef_plot.hurdle2.count[3,]), c(12,11,9:1) - .1, lwd = 2)
points(coef_plot.hurdle2.count[1,], c(12,11,9:1) - .1, pch = 24, bg = "white")
abline(v=0, lty = 2)
dev.off()



# hurdle models for figure 6 in the supplementary materials ---------------

hurdle3 <- hurdle(mentioned_negative ~ cjeu_pos_auto.numerical + observation_referring_ms + pos_auto.numerical + 
                    oral_observation + z.subject_matter.count + panel_size +
                    legact_primary + legact_secondary + z.gdp + largest_ms +
                    year_proceeding.factor,
                  data = mentions_submission_level, dist = "negbin")

coef.count3 <- summary(hurdle3)$coefficients$count
coef.zero3 <- summary(hurdle3)$coefficients$zero

hurdle4 <- hurdle(mentioned_negative ~ observation_referring_ms + pos_auto.numerical + 
                    oral_observation + z.subject_matter.count + panel_size +
                    legact_primary + legact_secondary + z.gdp + largest_ms +
                    year_proceeding.factor,
                  data = mentions_submission_level, dist = "negbin")

coef.count4 <- summary(hurdle4)$coefficients$count
coef.zero4 <- summary(hurdle4)$coefficients$zero

coef_plot.hurdle3.count <- hurdle_point_estimates(coef.count3, c("cjeu_pos_auto.numerical","observation_referring_ms","pos_auto.numerical",
                                                                 "oral_observation","z.subject_matter.count","panel_size",
                                                                 "legact_primary","legact_secondary","z.gdp","largest_ms"))

coef_plot.hurdle3.zero <- hurdle_point_estimates(coef.zero3, c("cjeu_pos_auto.numerical","observation_referring_ms","pos_auto.numerical",
                                                               "oral_observation","z.subject_matter.count","panel_size",
                                                               "legact_primary","legact_secondary","z.gdp","largest_ms"))

coef_plot.hurdle4.count <- hurdle_point_estimates(coef.count4, c("observation_referring_ms","pos_auto.numerical",
                                                                 "oral_observation","z.subject_matter.count","panel_size",
                                                                 "legact_primary","legact_secondary","z.gdp","largest_ms"))

coef_plot.hurdle4.zero <- hurdle_point_estimates(coef.zero4, c("observation_referring_ms","pos_auto.numerical",
                                                               "oral_observation","z.subject_matter.count","panel_size",
                                                               "legact_primary","legact_secondary","z.gdp","largest_ms"))


# plot coefficients for model 9 and 10
pdf(file = "figures/hurdle_submissions.pdf", height = 9, width = 10)
par(mfrow=c(2,1))
# zero model
par(mar=c(5,18,3,3))
plot("", axes = FALSE,
     xlim = c(-1.75,1.75),
     ylim = c(0.75,10.25),
     xlab = "Logistic regression coefficient estimates with 95% confidence intervals",
     ylab = "", pch = 19, cex = 1.5,
     main = "Zero model: Negative reference (Models 9 and 10)", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-1.75,1.75,.25), cex.axis = 1.2)
axis(2,at = seq(10,1,-1), las = 1, tick = TRUE, cex.axis = 1.2,
     label = c("CJEU position",
               "Case originated in MS",
               "Member State position",
               "Oral hearing",
               "Number of subject matters",
               "Panel size",
               "Issue concerns primary EU law",
               "Issue concerns secondary EU law",
               "Member State GDP per capita",
               "Largest Member States"))
abline(h = seq(10,1,-1), lty = 2, lwd = .5, col = "grey")
abline(v=0, lty = 2)
# m1
points(coef_plot.hurdle3.zero[1,], 10:1 + .1, pch = 19)
segments(t(coef_plot.hurdle3.zero[2,]), 10:1 + .1,
         t(coef_plot.hurdle3.zero[3,]), 10:1 + .1, lwd = 2)
# m2
segments(t(coef_plot.hurdle4.zero[2,]), 9:1 - .1,
         t(coef_plot.hurdle4.zero[3,]), 9:1 - .1, lwd = 2)
points(coef_plot.hurdle4.zero[1,], 9:1 - .1, pch = 24, bg = "white")
# count model
par(mar=c(5,18,3,3))
plot("", axes = FALSE,
     xlim = c(-1.25,1.25),
     ylim = c(0.75,10.25),
     xlab = "Negative binomial regression coefficient estimates with 95% confidence intervals",
     ylab = "", pch = 19, cex = 1.5,
     main = "Count model: Negative references (Models 9 and 10)", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-1.25,1.25,.25), cex.axis = 1.2)
axis(2,at = seq(10,1,-1), las = 1, tick = TRUE, cex.axis = 1.2,
     label = c("CJEU position",
               "Case originated in MS",
               "Member State position",
               "Oral hearing",
               "Number of subject matters",
               "Panel size",
               "Issue concerns primary EU law",
               "Issue concerns secondary EU law",
               "Member State GDP per capita",
               "Largest Member States"))
abline(h = seq(10,1,-1), lty = 2, lwd = .5, col = "grey")
# m1
points(coef_plot.hurdle3.count[1,], 10:1 + .1, pch = 19)
segments(t(coef_plot.hurdle3.count[2,]), 10:1 + .1,
         t(coef_plot.hurdle3.count[3,]), 10:1 + .1, lwd = 2)
# m2
segments(t(coef_plot.hurdle4.count[2,]), 9:1 - .1,
         t(coef_plot.hurdle4.count[3,]), 9:1 - .1, lwd = 2)
points(coef_plot.hurdle4.count[1,], 9:1 - .1, pch = 24, bg = "white")
abline(v=0, lty = 2)
dev.off()



# multinomial regression for figure 7 in the supplementary materia --------

# multinomial regression with cjeu position
ml_issue_cjeu <- multinom(mentioned_cat ~ cjeu_pos_auto.numerical + z.ms_pos_auto0 + z.ms_pos_auto1 + ag_cjeu_disagreement + ag_com_disagreement +
                            z.sources_case_law.ag + z.sources_primary.ag + z.sources_secondary.ag +
                            z.duration_since_ag + panel_size +
                            number_of_observations + legact_primary + legact_secondary +
                            administrative + agriculture + criminal + customs + environment + healthcare + intellectual_property +
                            labour + marketing + migration + procedure + procurement + no_area_of_law + social + tax + transport +
                            z.subject_matter.count + year_proceeding.factor, data = mentions_issue_level)
ml_issue_cjeu_coefs <- summary(ml_issue_cjeu)$coefficients[,2:14]
ml_issue_cjeu_coefs[1,]
ml_issue_cjeu_ses <- summary(ml_issue_cjeu)$standard.errors[,2:14]

ml_issue_cjeu_hi <- ml_issue_cjeu_coefs + 1.96 * ml_issue_cjeu_ses
ml_issue_cjeu_lo <- ml_issue_cjeu_coefs - 1.96 * ml_issue_cjeu_ses

# multinomial regression w/o cjeu position
ml_issue_no_cjeu <- multinom(mentioned_cat ~ z.ms_pos_auto0 + z.ms_pos_auto1 + ag_com_disagreement +
                               z.sources_case_law.ag + z.sources_primary.ag + z.sources_secondary.ag +
                               z.duration_since_ag + panel_size +
                               number_of_observations + legact_primary + legact_secondary +
                               administrative + agriculture + criminal + customs + environment + healthcare + intellectual_property +
                               labour + marketing + migration + procedure + procurement + no_area_of_law + social + tax + transport +
                               z.subject_matter.count + year_proceeding.factor, data = mentions_issue_level)
ml_issue_no_cjeu_coefs <- summary(ml_issue_no_cjeu)$coefficients[,2:12]
ml_issue_no_cjeu_ses <- summary(ml_issue_no_cjeu)$standard.errors[,2:12]

ml_issue_no_cjeu_hi <- ml_issue_no_cjeu_coefs + 1.96 * ml_issue_no_cjeu_ses
ml_issue_no_cjeu_lo <- ml_issue_no_cjeu_coefs - 1.96 * ml_issue_no_cjeu_ses

# plot coefficient estimates for each outcome category

ylabels <- c("CJEU position",
             "MS favouring no restrictions",
             "MS favouring restrictions",
             "Conflicting positions: AG/CJEU",
             "Conflicting positions: AG/Commission",
             "Sources cited by AG: Case law",
             "Sources cited by AG: Primary law",
             "Sources cited by AG: Secondary law",
             "Deliberation time",
             "Panel size",
             "Number of observations",
             "Issue concerns primary EU law",
             "Issue concerns secondary EU law")


pdf(file = "figures/multinomial_issue_top.pdf", height = 6, width = 10)
par(mfrow = c(1,2))
#--- negative
par(mar=c(2,2,2,0.5), oma = c(2, 13, 0.2, 0.2))
plot("", axes = FALSE,
     xlim = c(-1.75,1.75),
     ylim = c(0.75,13.25),
     xlab = "",
     ylab = "", pch = 19, cex = 1.5,
     main = "Negative reference", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-1.75,1.75,.25), cex.axis = 1.2)
axis(2,at = seq(13,1,-1), las = 1, tick = TRUE, cex.axis = 1.2,
     label = rep("",13))
abline(h = seq(13,1,-1), lty = 2, lwd = .5, col = "grey")
# with cjeu
points(ml_issue_cjeu_coefs[1,], 13:1 + 0.1, pch = 19)
segments(ml_issue_cjeu_hi[1,], 13:1 + 0.1,
         ml_issue_cjeu_lo[1,], 13:1 + 0.1, lwd = 2)
# w/o cjeu
segments(ml_issue_no_cjeu_hi[1,], c(12,11,9:1) - 0.1,
         ml_issue_no_cjeu_lo[1,], c(12,11,9:1) - 0.1, lwd = 2)
points(ml_issue_no_cjeu_coefs[1,], c(12,11,9:1) - 0.1, pch = 24, bg = "white")
abline(v=0, lty = 2)
for (i in length(ylabels):1){
  text(-3.8, i, rev(ylabels)[i], xpd=NA)
}

#--- neutral
par(mar=c(2,2,2,0.5))
plot("", axes = FALSE,
     xlim = c(-1.75,1.75),
     ylim = c(0.75,13.25),
     xlab = "",
     ylab = "", pch = 19, cex = 1.5,
     main = "Neutral reference", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-1.75,1.75,.25), cex.axis = 1.2)
abline(h = seq(13,1,-1), lty = 2, lwd = .5, col = "grey")
# with cjeu
points(ml_issue_cjeu_coefs[2,], 13:1 + 0.1, pch = 19)
segments(ml_issue_cjeu_hi[2,], 13:1 + 0.1,
         ml_issue_cjeu_lo[2,], 13:1 + 0.1, lwd = 2)
# w/o cjeu
segments(ml_issue_no_cjeu_hi[2,], c(12,11,9:1) - 0.1,
         ml_issue_no_cjeu_lo[2,], c(12,11,9:1) - 0.1, lwd = 2)
points(ml_issue_no_cjeu_coefs[2,], c(12,11,9:1) - 0.1, pch = 24, bg = "white")
abline(v=0, lty = 2)
dev.off()


pdf(file = "figures/multinomial_issue_bottom.pdf", height = 6, width = 10)
par(mfrow = c(1,2))
#--- positive
par(mar=c(2,2,2,0.5), oma = c(2, 13, 0.2, 0.2))
plot("", axes = FALSE,
     xlim = c(-1.75,1.75),
     ylim = c(0.75,13.25),
     xlab = "",
     ylab = "", pch = 19, cex = 1.5,
     main = "Positive reference", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-1.75,1.75,.25), cex.axis = 1.2)
abline(h = seq(13,1,-1), lty = 2, lwd = .5, col = "grey")
# with cjeu
points(ml_issue_cjeu_coefs[3,], 13:1 + 0.1, pch = 19)
segments(ml_issue_cjeu_hi[3,], 13:1 + 0.1,
         ml_issue_cjeu_lo[3,], 13:1 + 0.1, lwd = 2)
# w/o cjeu
segments(ml_issue_no_cjeu_hi[3,], c(12,11,9:1) - 0.1,
         ml_issue_no_cjeu_lo[3,], c(12,11,9:1) - 0.1, lwd = 2)
points(ml_issue_no_cjeu_coefs[3,], c(12,11,9:1) - 0.1, pch = 24, bg = "white")
abline(v=0, lty = 2)
for (i in length(ylabels):1){
  text(-3.8, i, rev(ylabels)[i], xpd=NA)
}

#--- positive & negative
par(mar=c(2,2,2,0.5))
plot("", axes = FALSE,
     xlim = c(-1.75,1.75),
     ylim = c(0.75,13.25),
     xlab = "",
     ylab = "", pch = 19, cex = 1.5,
     main = "Positive and negative reference", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-1.75,1.75,.25), cex.axis = 1.2)
abline(h = seq(13,1,-1), lty = 2, lwd = .5, col = "grey")
# with cjeu
points(ml_issue_cjeu_coefs[4,], 13:1 + 0.1, pch = 19)
segments(ml_issue_cjeu_hi[4,], 13:1 + 0.1,
         ml_issue_cjeu_lo[4,], 13:1 + 0.1, lwd = 2)
# w/o cjeu
segments(ml_issue_no_cjeu_hi[4,], c(12,11,9:1) - 0.1,
         ml_issue_no_cjeu_lo[4,], c(12,11,9:1) - 0.1, lwd = 2)
points(ml_issue_no_cjeu_coefs[4,], c(12,11,9:1) - 0.1, pch = 24, bg = "white")
abline(v=0, lty = 2)
dev.off()



# multinomial regression for figure 8 in the supplementary materia --------

# multinomial regression with cjeu position
ml_submission_cjeu <- multinom(mentioned_cat ~ cjeu_pos_auto.numerical + observation_referring_ms + pos_auto.numerical + 
                                 oral_observation + z.subject_matter.count + panel_size +
                                 legact_primary + legact_secondary + z.gdp + largest_ms, data = mentions_submission_level)
ml_submission_cjeu_coefs <- summary(ml_submission_cjeu)$coefficients[,2:11]
ml_submission_cjeu_ses <- summary(ml_submission_cjeu)$standard.errors[,2:11]

ml_submission_cjeu_hi <- ml_submission_cjeu_coefs + 1.96 * ml_submission_cjeu_ses
ml_submission_cjeu_lo <- ml_submission_cjeu_coefs - 1.96 * ml_submission_cjeu_ses

# multinomial regression w/o cjeu position
ml_submission_no_cjeu <- multinom(mentioned_cat ~ observation_referring_ms + pos_auto.numerical + 
                                    oral_observation + z.subject_matter.count + panel_size +
                                    legact_primary + legact_secondary + z.gdp + largest_ms, data = mentions_submission_level)
ml_submission_no_cjeu_coefs <- summary(ml_submission_no_cjeu)$coefficients[,2:10]
ml_submission_no_cjeu_ses <- summary(ml_submission_no_cjeu)$standard.errors[,2:10]

ml_submission_no_cjeu_hi <- ml_submission_no_cjeu_coefs + 1.96 * ml_submission_no_cjeu_ses
ml_submission_no_cjeu_lo <- ml_submission_no_cjeu_coefs - 1.96 * ml_submission_no_cjeu_ses


ylabels <- c("CJEU position",
             "Case originated in MS",
             "Member State position",
             "Oral hearing",
             "Number of subject matters",
             "Panel size",
             "Issue concerns primary EU law",
             "Issue concerns secondary EU law",
             "Member State GDP per capita",
             "Largest Member States")


pdf(file = "figures/multinomial_submission_top.pdf", height = 6, width = 10)
par(mfrow = c(1,2))
#--- negative
par(mar=c(2,2,2,0.5), oma = c(2, 13, 0.2, 0.2))
plot("", axes = FALSE,
     xlim = c(-1.75,1.75),
     ylim = c(0.75,10.25),
     xlab = "",
     ylab = "", pch = 19, cex = 1.5,
     main = "Negative reference", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-1.75,1.75,.25), cex.axis = 1.2)
axis(2,at = seq(10,1,-1), las = 1, tick = TRUE, cex.axis = 1.2,
     label = rep("",10))
abline(h = seq(10,1,-1), lty = 2, lwd = .5, col = "grey")
# with cjeu
points(ml_submission_cjeu_coefs[1,], 10:1 + 0.1, pch = 19)
segments(ml_submission_cjeu_hi[1,], 10:1 + 0.1,
         ml_submission_cjeu_lo[1,], 10:1 + 0.1, lwd = 2)
# w/o cjeu
segments(ml_submission_no_cjeu_hi[1,], 9:1 - 0.1,
         ml_submission_no_cjeu_lo[1,], 9:1 - 0.1, lwd = 2)
points(ml_submission_no_cjeu_coefs[1,], 9:1 - 0.1, pch = 24, bg = "white")
abline(v=0, lty = 2)
for (i in length(ylabels):1){
  text(-3.8, i, rev(ylabels)[i], xpd=NA)
}

#--- neutral
par(mar=c(2,2,2,0.5))
plot("", axes = FALSE,
     xlim = c(-1.75,1.75),
     ylim = c(0.75,10.25),
     xlab = "",
     ylab = "", pch = 19, cex = 1.5,
     main = "Neutral reference", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-1.75,1.75,.25), cex.axis = 1.2)
abline(h = seq(10,1,-1), lty = 2, lwd = .5, col = "grey")
# with cjeu
points(ml_submission_cjeu_coefs[2,], 10:1 + 0.1, pch = 19)
segments(ml_submission_cjeu_hi[2,], 10:1 + 0.1,
         ml_submission_cjeu_lo[2,], 10:1 + 0.1, lwd = 2)
# w/o cjeu
segments(ml_submission_no_cjeu_hi[2,], 9:1 - 0.1,
         ml_submission_no_cjeu_lo[2,], 9:1 - 0.1, lwd = 2)
points(ml_submission_no_cjeu_coefs[2,], 9:1 - 0.1, pch = 24, bg = "white")
abline(v=0, lty = 2)
dev.off()


pdf(file = "figures/multinomial_submission_bottom.pdf", height = 6, width = 10)
par(mfrow = c(1,2))
#--- positive
par(mar=c(2,2,2,0.5), oma = c(2, 13, 0.2, 0.2))
plot("", axes = FALSE,
     xlim = c(-1.75,1.75),
     ylim = c(0.75,10.25),
     xlab = "",
     ylab = "", pch = 19, cex = 1.5,
     main = "Positive reference", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-1.75,1.75,.25), cex.axis = 1.2)
abline(h = seq(10,1,-1), lty = 2, lwd = .5, col = "grey")
# with cjeu
points(ml_submission_cjeu_coefs[3,], 10:1 + 0.1, pch = 19)
segments(ml_submission_cjeu_hi[3,], 10:1 + 0.1,
         ml_submission_cjeu_lo[3,], 10:1 + 0.1, lwd = 2)
# w/o cjeu
segments(ml_submission_no_cjeu_hi[3,], 9:1 - 0.1,
         ml_submission_no_cjeu_lo[3,], 9:1 - 0.1, lwd = 2)
points(ml_submission_no_cjeu_coefs[3,], 9:1 - 0.1, pch = 24, bg = "white")
abline(v=0, lty = 2)
for (i in length(ylabels):1){
  text(-3.8, i, rev(ylabels)[i], xpd=NA)
}

#--- positive & negative
par(mar=c(2,2,2,0.5))
plot("", axes = FALSE,
     xlim = c(-1.75,3),
     ylim = c(0.75,10.25),
     xlab = "",
     ylab = "", pch = 19, cex = 1.5,
     main = "Positive and negative reference", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-1.75,3,.5), cex.axis = 1.2)
abline(h = seq(10,1,-1), lty = 2, lwd = .5, col = "grey")
# with cjeu
points(ml_submission_cjeu_coefs[4,], 10:1 + 0.1, pch = 19)
segments(ml_submission_cjeu_hi[4,], 10:1 + 0.1,
         ml_submission_cjeu_lo[4,], 10:1 + 0.1, lwd = 2)
# w/o cjeu
segments(ml_submission_no_cjeu_hi[4,], 9:1 - 0.1,
         ml_submission_no_cjeu_lo[4,], 9:1 - 0.1, lwd = 2)
points(ml_submission_no_cjeu_coefs[4,], 9:1 - 0.1, pch = 24, bg = "white")
abline(v=0, lty = 2)
dev.off()
